% Calculate and plot similarity matrices
% Plot basis functions and energy in modes
clear all; close all; clc;

%species = {'Ciona'; 'Arbacia'; 'Lpictus'; 'Lvariegatus';...
%    'Spurpuratus'; 'Zebrafish'; 'Abalone'; 'Bull'; ...
%    'Human-low'; 'Human-high';};
%species = {'Abalone';'Lpictus';'Lvariegatus';'Spurpuratus'; ...
%    'Arbacia';'Ciona';'Zebrafish';'Bull';'Human-low';'Human-high';};

spec = 'Lpictus';

switch spec
    case 'Spurpuratus'
        species = {'Spurpuratus','SpurpuratusD1','SpurpuratusD2','SpurpuratusD3'};
        date = {'200617B'; '200617'; '200617';'200617'};
        trlens = [30322, 51811, 46756];
    case 'Ciona'
        species = {'Ciona','CionaD1','CionaD2'};
        date = {'200617B'; '200617'; '200617'};
        trlens = [34703, 23341];
    case 'Zebrafish'
        species = {'Zebrafish','ZebrafishD1','ZebrafishD2','ZebrafishD3'};
        date = {'200617B'; '200617'; '200617';'200617'};
        trlens = [4410, 11998, 3256];
    case 'Lpictus'
        species = {'Lpictus','LpictusD1','LpictusD2','LpictusD3','LpictusD4'};
        date = {'200617B'; '200617'; '200617';'200617';'200617'};
        trlens = [86921, 46617, 70071, 61486];
    case 'Lvariegatus'
        species = {'Lvariegatus','LvariegatusD1','LvariegatusD2','LvariegatusD3','LvariegatusD4'};
        date = {'200617B'; '200617'; '200617';'200617';'200617'};
        trlens = [36349, 11056, 18366, 35086];
end

use_optimization = 1; % 1 - use optimization, 0 - don't use optimization
dat = struct([]);
recalcAll = false;

savedate = ['200622-KD-Donorcomp' spec];

actualSpecies{1} = [spec '-All'];
for i=length(date):-1:2
actualSpecies{i} = [spec '-' num2str(i-1)];
end

% actualSpecies{1} = 'H. rufescens';
% actualSpecies{2} = 'L. pictus';
% actualSpecies{3} = 'L. variegatus';
% actualSpecies{4} = 'S. purpuratus';
% actualSpecies{5} = 'A. punctulata';
% actualSpecies{6} = 'C. intestinalis';
% actualSpecies{7} = 'D. rerio';
% actualSpecies{8} = 'B. taurus';
% actualSpecies{9} = 'H. sapiens';
% actualSpecies{10} = 'H. sapiens (viscous)';

symb = 'ovs^oooovp';

for x = 1 : length(species)
    
    filename = strcat(species{x},['_' date{x} '.mat']);
    load(filename)
    dat(x).Uall = U_temp;
    dat(x).Sall = S_temp;
    %if exist('procdata')
    %    dat(x).U = double(procdata(1).U);
    %    dat(x).S = double(procdata(1).S);
    %    clear procdata allruns U_temp thresh kfolds
    %else
    U_temp = cat(3,U_temp{:});
    
    stdU{x} = std(U_temp,[],3);
    U_temp = mean(U_temp,3);
    sEnd = size(S_temp{1},1);
    for i =1:length(S_temp)
        try
            S_temp{i} = S_temp{i}(1:sEnd,1:sEnd);
        catch
            sz = size(S_temp{i});
            S_temp{i} = [S_temp{i}(1:sEnd,1:sz(2)), zeros(sEnd,sEnd-sz(2))];
        end
        svi = diag(S_temp{i});
        dat(x).Sall{i} = double(cumsum(svi.^2) ./ sum(svi.^2));
    end
    S_temp = cat(3,S_temp{:});
    S_temp = mean(S_temp,3);
    sv = diag(S_temp);
    dat(x).U = double(U_temp);
    
    dat(x).S = double(cumsum(sv.^2) ./ sum(sv.^2));
    
    
    clear A_overthresh allruns U_temp thresh kfolds S_temp
    %end
end

%% Kinematic Error on all individual curves
%wMval = zeros(10,10,sum(1:49));
wMval = nan(length(species),length(species),100);

s = RandStream('mlfg6331_64');
for p = 1:100
    rands100(p,:) = datasample(s,1:50,2,'Replace',false);
end
%rands100 = randi(50,100,2);

for i = length(dat):-1:1
    for j = length(dat):-1:1
        counter = 0;
%         try
%             loadfile = load('200617-EL2-KD-Final_ALLwMerror.mat');
%             wMval(~isnan(loadfile.wMval)) = loadfile.wMval(~isnan(loadfile.wMval));
%         end
        try
            loadfile = load([savedate '_ALLwMerror.mat']);
            wMval(~isnan(loadfile.wMval)) = loadfile.wMval(~isnan(loadfile.wMval));
        end
        %for ind1 = 1:50
        tic
        for p=1:50
            try
                ind2 = rands100(p,1);
                ind1 = rands100(p,2);
                %for ind2 = (ind1+1):50
                counter = counter+1;
                U_{1} = dat(i).Uall{ind1}; U_{2} = dat(j).Uall{ind2};
                %S_{1} = dat(i).Sall{ind1}; S_{2} = dat(j).Sall{ind2};
                S_{1} = dat(i).Sall{ind1}; S_{2} = dat(j).Sall{ind2};
                if isnan(wMval(i,j,counter)) || recalcAll
                    wMval(i,j,counter) = calculateWMval(U_, S_);
                end
                
                while wMval(i,j,counter) > 1 || wMval(i,j,counter) < 0
                    [ind1, ind2] = datasample(s,1:50,2,'Replace',false);
                    U_{1} = dat(i).Uall{ind1}; U_{2} = dat(j).Uall{ind2};
                    %S_{1} = dat(i).Sall{ind1}; S_{2} = dat(j).Sall{ind2};
                    S_{1} = dat(i).Sall{ind1}; S_{2} = dat(j).Sall{ind2};
                    wMval(i,j,counter) = calculateWMval(U_, S_);
                end
                %end
            end
            tempvals = wMval(i,j,:);
            disp(['Mean through ' num2str(counter) ' of ' num2str(100) ': ' num2str(nanmedian(tempvals(tempvals~=0)))]);
            disp(['Std through ' num2str(counter) ' of ' num2str(100) ': ' num2str(nanstd(tempvals(tempvals~=0)))]);
            
            if mod(p,10)==0
                save([savedate '_ALLwMerror.mat'],'wMval');
                figure(1); imagesc(nanmedian(wMval,3)); axis image; colormap hot; colorbar;
                figure(2); imagesc(100-sum(isnan(wMval),3)); axis image; colormap gray; colorbar;
            end
            
        end
        toc
        
    end
end

%medWM = nanmedian(wMval,3);
% KD_ = medWM; 
% KD_(2,:) = medWM(4,:); KD_(:,2) = medWM(:,4);
% KD_(4,:) = medWM(2,:); KD_(:,4) = medWM(:,2);
% KD_(2,4) = medWM(4,2); KD_(4,2) = medWM(2,4);
% KD_(2,2) = medWM(4,4); KD_(4,4) = medWM(2,2);
% medWM = KD_;


%%
clear all; close all;

spec = 'Spurpuratus';
load(['200622-KD-Donorcomp' spec '_ALLwMerror.mat']);

switch spec
    case 'Spurpuratus'
        species = {'Spurpuratus','SpurpuratusD1','SpurpuratusD2','SpurpuratusD3'};
        date = {'200617'; '200617'; '200617';'200617'};
        trlens = [30322, 51811, 46756];
    case 'Ciona'
        species = {'Ciona','CionaD1','CionaD2'};
        date = {'200617'; '200617'; '200617'};
        trlens = [34703, 23341];
    case 'Zebrafish'
        species = {'Zebrafish','ZebrafishD1','ZebrafishD2','ZebrafishD3'};
        date = {'200617'; '200617'; '200617';'200617'};
        trlens = [4410, 11998, 3256];
    case 'Lpictus'
        species = {'Lpictus','LpictusD1','LpictusD2','LpictusD3','LpictusD4'};
        date = {'200617'; '200617'; '200617';'200617';'200617'};
        trlens = [86921, 46617, 70071, 61486];
    case 'Lvariegatus'
        species = {'Lvariegatus','LvariegatusD1','LvariegatusD2','LvariegatusD3','LvariegatusD4'};
        date = {'200617'; '200617'; '200617';'200617';'200617'};
        trlens = [36349, 11056, 18366, 35086];
end
actualSpecies{1} = [spec '-All'];
for i=length(date):-1:2
actualSpecies{i} = [spec '-' num2str(i-1)];
end


medWM = nanmedian(wMval,3);
wMij = mean(cat(3,medWM,medWM'),3);
figure; imagesc(wMij); axis image; colormap hot; colorbar;
set(gca,'xtick',[1:10],'xticklabel',actualSpecies,'ytick',[1:10],'yticklabel',actualSpecies);
xtickangle(90);
caxis([0 0.071]);

%figure; imagesc(wMij(1:9,1:9)); axis image; colormap hot; colorbar;

for i=1:size(wMij,2)
    for k=1:size(wMij,2)
    bpWMt{i}{k} = squeeze(wMval(i,k,:)); 
    [n(i,k),h(i,k)] = ranksum(squeeze(wMval(1,i,1:50)),squeeze(wMval(1,k,1:50)));
    end
end
% ranksum(squeeze(wMval(1,2,1:50)),squeeze(wMval(1,3,1:50)))
figure; violin(bpWMt{1},[],'labels',actualSpecies,'Color',[0.0 0.5 0])
figure; subplot(1,2,1); imagesc(n); axis image; subplot(1,2,2); imagesc(h); axis image;


%mean donor-donor distance weighted by track
numdonors = size(wMval,2)-1;
dg = eye(numdonors)==1;
meanDD = 0;
temptrlens = 0;
for i=1:numdonors
    for j=1:numdonors
       if i~=j
           temptrlens = temptrlens + (trlens(i)+trlens(j))/2;
           meanDD = meanDD + wMij(i+1,j+1)*(trlens(i)+trlens(j))/2;
       end
    end
end
meanDD = meanDD/temptrlens;

meanDS = 0;
for j=1:numdonors
       meanDS = meanDS + wMij(1,j+1)*trlens(j);
end
meanDS = meanDS/sum(trlens);




